#!/bin/bash -e

function info() {
echo Usage: `basename $0` snp.dbsnp.af.vcf
exit 65
}

while getopts  ":db:p:" opts
do
        case  $opts  in
		b) asemmbly=$OPTARG;;
		d) anno_dbsnp=true;;
		p) out_prefix=$OPTARG;;
		*) info;;
        esac
done
shift $(($OPTIND - 1))

test -z "$1" && info


. /mnt/ilustre/app/medical/tools/.var

sample_name=$(grep -m1 '^#C' $1|awk '{print $NF}')

af_min=0.01
af_min_percent=`echo "$af_min * 100"|bc`
adp_min=6
adp_min_alt=2
mq_min=20.0



echo
echo
echo gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o $out_prefix.snp.filtered.vcf \
	--variant $1 \
	--filterExpression "FS > 60.0" \
	--filterName "FS_more_than_60" \
	--filterExpression "MQ < $mq_min" \
	--filterName "MQ_less_than_$mq_min" \
	--filterExpression "HaplotypeScore > 200.0" \
	--filterName "HaplotypeScore_more_than_200" \
	--filterExpression "af < $af_min" \
	--filterName "Allele_frequency_less_than_$af_min_percent%" \
	--filterExpression "DP < $adp_min" \
	--filterName "Depth_less_than_$adp_min"



slt="vc.getGenotype(\"$sample_name\").getAD().1 > $adp_min_alt && vc.isNotFiltered()"


echo
echo
echo gatk SelectVariants
echo $slt
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T SelectVariants \
	--variant $out_prefix.snp.filtered.vcf \
	-o $out_prefix.filtered.pass.vcf \
	-select "$slt"


cat $out_prefix.filtered.pass.vcf|grep -v "#" | awk '{if($1~/^[1-9]$/){print 0$1$2$3"/"$NF}else{print $1$2$3"/"$NF}}' | cut -d':' -f1|awk -F "/" '{if($2==$3){print $1 > "'"$out_prefix"'.hom.txt"}else{print $1 > "'"$out_prefix"'.het.txt"}}'

cat $out_prefix.hom.txt|grep rs > $out_prefix.hom.rs.txt
cat $out_prefix.het.txt|grep rs > $out_prefix.het.rs.txt
